home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / clustalw.src / calcgapcoeff.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  8.3 KB  |  340 lines  |  [TEXT/R*ch]
  1. #include <stdio.h>
  2. #include <ctype.h>
  3. #include <stdlib.h>
  4. #include <string.h>
  5. #include "clustalw.h"
  6.  
  7. /*
  8.  *   Prototypes
  9.  */
  10. extern void *ckalloc(size_t);
  11. extern void ckfree(void *);
  12.  
  13. void calc_p_penalties(char **aln, int n, int fs, int ls, int *weight);
  14. void calc_h_penalties(char **aln, int n, int fs, int ls, int *weight);
  15. int  local_penalty(int penalty, int n, int *pweight, int *hweight);
  16. void calc_gap_coeff(char **alignment, int *gaps, int **profile,
  17.         int first_seq, int last_seq, int prf_length, int gapcoef, int lencoef);
  18.  
  19. /*
  20.  *   Global variables
  21.  */
  22.  
  23. extern int gap_dist;
  24. extern int max_aa;
  25. extern int debug;
  26. extern Boolean dnaflag;
  27. extern Boolean use_endgaps;
  28. extern Boolean no_hyd_penalties, no_pref_penalties;
  29. extern char hyd_residues[];
  30. extern char *amino_acid_codes;
  31.  
  32. char   pr[] =     {'A' , 'C', 'D', 'E', 'F', 'G', 'H', 'K', 'I', 'L',
  33.                    'M' , 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'Y', 'W'};
  34. int    pas_op[] = { 87, 87,104, 69, 80,139,100,104, 68, 79,
  35.                     71,137,126, 93,128,124,111, 75,100, 77};
  36. int    pas_op2[] ={ 88, 57,111, 98, 75,126, 95, 97, 70, 90,
  37.                     60,122,110,107, 91,125,124, 81,106, 88};
  38. int    pal_op[] = { 84, 69,128, 78, 88,176, 53, 95, 55, 49,
  39.                     52,148,147,100, 91,129,105, 51,128, 88};
  40.  
  41. float reduced_gap = 0.3;
  42.  
  43. void calc_gap_coeff(char **alignment, int *gaps, int **profile,
  44.          int first_seq, int last_seq, int prf_length, int gapcoef, int lencoef)
  45. {
  46.  
  47.    char c;
  48.    int i, j;
  49.    int is, ie;
  50.    static int numseq;
  51.    static int *gap_pos;
  52.    static int *p_weight, *h_weight;
  53.    static float scale;
  54.  
  55.    numseq = last_seq - first_seq;
  56.  
  57.    for (j=0; j<prf_length; j++)
  58.         gaps[j] = 0;
  59.           
  60.    for (i=first_seq; i<last_seq; i++)
  61.      {
  62. /*
  63.    Include end gaps as gaps ?
  64. */
  65.         is = 0;
  66.         ie = prf_length;
  67.         if (use_endgaps == FALSE)
  68.         {
  69.           for (j=0; j<prf_length; j++)
  70.             {
  71.               c = alignment[i][j];
  72.               if ((c < 0) || (c > max_aa))
  73.                  is++;
  74.               else
  75.                  break;
  76.             }
  77.           for (j=prf_length-1; j>=0; j--)
  78.             {
  79.               c = alignment[i][j];
  80.               if ((c < 0) || (c > max_aa))
  81.                  ie--;
  82.               else
  83.                  break;
  84.             }
  85.         }
  86.  
  87.         for (j=is; j<ie; j++)
  88.           {
  89.               if ((alignment[i][j] < 0) || (alignment[i][j] > max_aa))
  90.                  gaps[j]++;
  91.           }
  92.      }
  93.  
  94.    if ((!dnaflag) && (no_pref_penalties == FALSE))
  95.      {
  96.         p_weight = (int *) ckalloc( (prf_length+2) * sizeof (int) );
  97.         calc_p_penalties(alignment, prf_length, first_seq, last_seq, p_weight);
  98.      }
  99.  
  100.    if ((!dnaflag) && (no_hyd_penalties == FALSE))
  101.      {
  102.         h_weight = (int *) ckalloc( (prf_length+2) * sizeof (int) );
  103.         calc_h_penalties(alignment, prf_length, first_seq, last_seq, h_weight);
  104.      }
  105.  
  106.    gap_pos = (int *) ckalloc( (prf_length+2) * sizeof (int) );
  107. /*
  108.     mark the residues close to an existing gap (set gaps[i] = -ve)
  109. */
  110.    if (dnaflag || (gap_dist <= 0))
  111.      {
  112.        for (i=0;i<prf_length;i++) gap_pos[i] = gaps[i];
  113.      }
  114.    else
  115.      {
  116.        i=0;
  117.        while (i<prf_length)
  118.          {
  119.             if (gaps[i] <= 0)
  120.               {
  121.                  gap_pos[i] = gaps[i];
  122.                  i++;
  123.               }
  124.             else 
  125.               {
  126.                  for (j = -gap_dist+1; j<0; j++)
  127.                   {
  128.                    if ((i+j>=0) && (i+j<prf_length) &&
  129.                        ((gaps[i+j] == 0) || (gaps[i+j] < j))) gap_pos[i+j] = j;
  130.                   }
  131.                  while (gaps[i] > 0)
  132.                     {
  133.                        if (i>=prf_length) break;
  134.                        gap_pos[i] = gaps[i];
  135.                        i++;
  136.                     }
  137.                  for (j = 0; j<gap_dist; j++)
  138.                   {
  139.                    if (gaps[i+j] > 0) break;
  140.                    if ((i+j>=0) && (i+j<prf_length) && 
  141.                        ((gaps[i+j] == 0) || (gaps[i+j] < -j))) gap_pos[i+j] = -j-1;
  142.                   }
  143.                  i += j;
  144.               }
  145.          }
  146.      }
  147.  
  148. if (debug>1)
  149. {
  150. fprintf(stderr,"gap open %d gap ext %d\n",(pint)gapcoef,(pint)lencoef);
  151.   for(i=0;i<prf_length;i++) fprintf(stderr,"%d ", (pint)gaps[i]);
  152.   fprintf(stderr,"\n");
  153.   for(i=0;i<prf_length;i++) fprintf(stderr,"%d ", (pint)gap_pos[i]);
  154.   fprintf(stderr,"\n");
  155. }
  156.  
  157.    for (j=0;j<prf_length; j++)
  158.      {
  159.         if (gap_pos[j] <= 0)
  160.           {
  161. /*
  162.     apply residue-specific and hydrophilic gap penalties.
  163. */
  164.          if (!dnaflag) {
  165.                   profile[j+1][GAPCOL] = local_penalty(gapcoef, j,
  166.                                                    p_weight, h_weight);
  167.                   profile[j+1][LENCOL] = lencoef;
  168.          }
  169.          else {
  170.                   profile[j+1][GAPCOL] = gapcoef;
  171.                   profile[j+1][LENCOL] = lencoef;
  172.          }
  173.  
  174. /*
  175.     increase gap penalty near to existing gaps.
  176. */
  177.              if (gap_pos[j] < 0)
  178.                 {
  179.                     profile[j+1][GAPCOL] *= 2.0+2.0*(gap_dist+gap_pos[j])/gap_dist;
  180.                 }
  181.  
  182.  
  183.           }
  184.         else
  185.           {
  186.              scale = ((float)(numseq-gaps[j])/(float)numseq) * reduced_gap;
  187.              profile[j+1][GAPCOL] = scale*gapcoef;
  188.              profile[j+1][LENCOL] = 0.5 * lencoef;
  189.           }
  190.      }
  191.  
  192.    profile[0][GAPCOL] = 0;
  193.    profile[0][LENCOL] = 0;
  194.  
  195.    profile[prf_length][GAPCOL] = 0;
  196.    profile[prf_length][LENCOL] = 0;
  197.  
  198. if (debug>0)
  199. {
  200.   fprintf(stderr,"Opening penalties:\n");
  201.   for(i=0;i<prf_length;i++) fprintf(stderr,"%d ", (pint)profile[i+1][GAPCOL]);
  202.   fprintf(stderr,"\n");
  203. }
  204. if (debug>0)
  205. {
  206.   fprintf(stderr,"Extension penalties:\n");
  207.   for(i=0;i<prf_length;i++) fprintf(stderr,"%d ", (pint)profile[i+1][LENCOL]);
  208.   fprintf(stderr,"\n");
  209. }
  210.    if ((!dnaflag) && (no_pref_penalties == FALSE))
  211.         ckfree((void *)p_weight);
  212.  
  213.    if ((!dnaflag) && (no_hyd_penalties == FALSE))
  214.         ckfree((void *)h_weight);
  215.  
  216.    ckfree((void *)gap_pos);
  217. }              
  218.             
  219. void calc_p_penalties(char **aln, int n, int fs, int ls, int *weight)
  220. {
  221.  
  222.   int i,j, k, ix;
  223.   int numseq;
  224.  
  225.   numseq = ls - fs;
  226.   for (i=0;i<n;i++)
  227.     {
  228.       weight[i] = 0;
  229.       for (k=fs;k<ls;k++)
  230.         {
  231.            for (j=0;j<22;j++)
  232.              {
  233.                 ix = aln[k][i];
  234.                 if ((ix < 0) || (ix > max_aa)) continue;
  235.                 if (amino_acid_codes[ix] == pr[j])
  236.                   {
  237.                     weight[i] += (200-pas_op[j]);
  238. /*
  239.                     weight[i] += (200-(pas_op[j]+pas_op[j+1])/2);
  240.                     weight[i] += ((200-pal_op[j])-100)/2+100;
  241.                     weight[i] += (200-pal_op[j]);
  242. */
  243.                     break;
  244.                   }
  245.              }
  246.         }
  247.       weight[i] /= numseq;
  248.     }
  249.  
  250. }
  251.             
  252. void calc_h_penalties(char **aln, int n, int fs, int ls, int *weight)
  253. {
  254.  
  255. /*
  256.    weight[] is the length of the hydrophilic run of residues.
  257. */
  258.  
  259.   int nh, h;
  260.   int i,j, ix, k, e, s, *hyd;
  261.   float scale;
  262.  
  263.   hyd = (int *)ckalloc((n+2) * sizeof(int));
  264.   nh = strlen(hyd_residues);
  265.   for (i=0;i<n;i++)
  266.      weight[i] = 0;
  267.  
  268.   for (k=fs;k<ls;k++)
  269.     {
  270.        for (i=0;i<n;i++)
  271.          {
  272.              hyd[i] = 0;
  273.              for (j=0;j<nh;j++)
  274.                 {
  275.                    ix = aln[k][i];
  276.                    if ((ix < 0) || (ix > max_aa)) continue;
  277.                    if (amino_acid_codes[ix] == hyd_residues[j])
  278.                       {
  279.                          hyd[i] = 1;
  280.                          break;
  281.                       }
  282.                 }
  283.           }
  284.        i = 0;
  285.        while (i < n)
  286.          {
  287.             if (hyd[i] == 0) i++;
  288.             else
  289.               {
  290.                  s = i;
  291.                  while ((hyd[i] != 0) && (i<n)) i++;
  292.                  e = i;
  293.                  if (e-s > 3)
  294.                     for (j=s; j<e; j++) weight[j] += 100;
  295.               }
  296.          }
  297.     }
  298.  
  299.   scale = ls - fs;
  300.   for (i=0;i<n;i++)
  301.      weight[i] /= scale;
  302.  
  303.   ckfree((void *)hyd);
  304.  
  305. if (debug>1)
  306. {
  307.   for(i=0;i<n;i++) fprintf(stderr,"%d ", (pint)weight[i]);
  308.   fprintf(stderr,"\n");
  309. }
  310.  
  311. }
  312.             
  313. int local_penalty(int penalty, int n, int *pweight, int *hweight)
  314. {
  315.  
  316.   int i,j, h = FALSE;
  317.   float gw;
  318.  
  319.   if (dnaflag) return(1);
  320.  
  321.   gw = 1.0;
  322.   if (no_hyd_penalties == FALSE)
  323.     {
  324.         if (hweight[n] > 0)
  325.          {
  326.            gw *= reduced_gap;
  327.            h = TRUE;
  328.          }
  329.     }
  330.   if ((no_pref_penalties == FALSE) && (h==FALSE))
  331.     {
  332.        gw *= ((float)pweight[n]/100.0);
  333.     }
  334.  
  335.   gw *= penalty;
  336.   return((int)gw);
  337.  
  338. }
  339.  
  340.